Update on assignments

• Assignment 2 due by 14 Oct

• Self-assessment 2 due 21 Oct

• Resubmits

• Final Project

Objectives

By the end of today you should be able to:

• Create new features based on topological relationships

• Use topological subsetting to reduce features

• Use spatial joins to add attributes based on location

What is spatial analysis?

“The process of examining the locations, attributes, and relationships of features in spatial data through overlay and other analytical techniques in order to address a question or gain useful knowledge. Spatial analysis extracts or creates new information from spatial data”.

Workflows for spatial analysis

Databases and Attributes

Attributes based on geometry and location (measures)

• Attributes like area and length can be useful for a number of analyses

  • Estimates of ‘effort’ in sampling designs
  • Offsets for modeling rates (e.g., Poisson regression)

• Need to assign the result of the function to a column in data frame (e.g., $, mutate, and summarize)

• Often useful to test before assigning

Estimating area

nz.sf <- nz %>% 
  mutate(area = st_area(nz))
head(nz.sf$area, 3)

nz.sf$areakm <- units::set_units(st_area(nz), km^2)
head(nz.sf$areakm, 3)

Estimating Density in Polygons

Estimating Density in Polygons

nz.df <- nz %>% 
mutate(counts = lengths(st_intersects(., random_nz)),
       area = st_area(nz),
       density = counts/area)
head(st_drop_geometry(nz.df[,7:10]))

Estimating Density in Polygons

Estimating Distance

• As a covariate

• For use in covariance matrices

• As a means of assigning connections in networks

Estimating Single Point Distance

Estimating Single Point Distance

• Subsetting y into a single feature

canterbury = nz %>% filter(Name == "Canterbury")
canterbury_height = nz_height[canterbury, ]
co = filter(nz, grepl("Canter|Otag", Name))
st_distance(nz_height[1:3, ], co)

Estimating Single Point Distance

• Using nearest neighbor distances

ua <- urban_areas(cb = FALSE, progress_bar = FALSE) %>% 
  filter(., UATYP10 == "U") %>% 
  filter(., str_detect(NAME10, "ID")) %>% 
  st_transform(., crs=2163)

#get index of nearest ID city
nearest <-  st_nearest_feature(ua)
#estimate distance
(dist = st_distance(ua, ua[nearest,], by_element=TRUE))

Topological Subsetting

ctby_height <-  nz_height[canterbury, ]

Topological Subsetting

st_intersects(nz_height, co, 
              sparse = TRUE)[1:3] 

lengths(st_intersects(nz_height, 
                      co, sparse = TRUE))[1:3] > 0

Topological Subsetting

• The sparse option controls how the results are returned
• We can then find out if one or more elements satisfies the criteria

Using sparse=FALSE

st_intersects(nz_height, co, sparse = FALSE)[1:3,] 

apply(st_intersects(nz_height, co, sparse = FALSE), 1,any)[1:3]

Topological Subsetting

canterbury_height3 = nz_height %>%
  filter(st_intersects(x = ., y = canterbury, sparse = FALSE))

Spatial Joins

• sf package provides st_join for vectors

• Allows joins based on the predicates (st_intersects, st_touches, st_within_distance, etc.)

• Default is a left join

Spatial Joins

set.seed(2018)
(bb = st_bbox(world)) # the world's bounds
#>   xmin   ymin   xmax   ymax 
#> -180.0  -89.9  180.0   83.6
random_df = data.frame(
  x = runif(n = 10, min = bb[1], max = bb[3]),
  y = runif(n = 10, min = bb[2], max = bb[4])
)
random_points = random_df |> 
  st_as_sf(coords = c("x", "y")) |> # set coordinates
  st_set_crs("EPSG:4326") # set geographic CRS

random_joined = st_join(random_points, world["name_long"])

Spatial Joins

• Sometimes we may want to be less restrictive

• Just because objects don’t touch doesn’t mean they don’t relate to each other

• Can use predicates in st_join

• Remember that default is left_join (so the number of records can grow if multiple matches)

Spatial Joins

any(st_touches(cycle_hire, cycle_hire_osm, sparse = FALSE))
z = st_join(cycle_hire, cycle_hire_osm, st_is_within_distance, dist = 20)
nrow(cycle_hire)
nrow(z)

Extending Joins

intersect_pct <- st_intersection(nc, tr_buff) %>% 
   mutate(intersect_area = st_area(.)) %>%   # create new column with shape area
   dplyr::select(NAME, intersect_area) %>%   # only select columns needed to merge
   st_drop_geometry()

nc <- mutate(nc, county_area = st_area(nc))

# Merge by county name
nc <- merge(nc, intersect_pct, by = "NAME", all.x = TRUE)

# Calculate coverage
nc <- nc %>% 
   mutate(coverage = as.numeric(intersect_area/county_area))

Extending Joins